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Abstract — We present a method for automated segmentation 
of the vasculature in retinal images. The method produces 
segmentations by classifying each image pixel as vessel or non- 
vessel, based on the pixel's feature vector. Feature vectors are 
composed of the pixel's intensity and continuous two-dimensional 
Morlet wavelet transform responses taken at multiple scales. The 
Morlet wavelet is capable of tuning to specific frequencies, thus 
allowing noise filtering and vessel enhancement in a single step. 
We use a Bayesian classifier with class-conditional probability 
density functions (likelihoods) described as Gaussian mixtures, 
yielding a fast classification, while being able to model complex 
decision surfaces and compare its performance with the linear 
minimum squared error classifier. The probability distributions 
are estimated based on a training set of labeled pixels ob- 
tained from manual segmentations. The method's performance 
is evaluated on publicly available DRIVE [1] and STARE [2] 
databases of manually labeled non-mydriatic images. On the 
DRIVE database, it achieves an area under the receiver operating 
characteristic (ROC) curve of 0.9598, being slightly superior than 
that presented by the method of Staal et al. [1]. 

Index Terms — Fundus, Morlet, pattern classification, retina, 
vessel segmentation, wavelet. 



I. Introduction 



OPTIC fundus (Fig. |l(a)| i assessment has been widely 
used by the medical community for diagnosing vas- 
cular and non-vascular pathology. Inspection of the retinal 
vasculature may reveal hypertension, diabetes, arteriosclerosis, 
cardiovascular disease and stroke [3]. Diabetic retinopathy is 
a major cause of adult blindness due to changes in blood 
vessel structure and distribution such as new vessel growth 
(proliferative diabetic retinopathy) and requires laborious anal- 
ysis from a specialist [4]. Endeavoring to reduce the effect 
of proliferative diabetic retinopathy includes obtaining and 
analyzing images of the optic fundus at regular intervals such 
as every six months to a year. Early recognition of changes 
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(a) Inverted green channel of 
non-mydriatic fundus image. 



(b) Pre-processed image with ex- 
tended border. The original im- 
age limit is indicated for illustra- 
tion. 



Fig. 1. Fundus image pre-processing for removing undesired border effects. 



to the blood vessel patterns can prevent major vision loss as 
early intervention becomes possible [5], [6]. 

To provide the opportunity for initial assessment to be 
carried out by community health workers, computer based 
analysis has been introduced, which includes assessment of 
the presence of microaneurysms and changes in the blood 
flow/vessel distribution due to either vessel narrowing, com- 
plete occlusions or new vessel growth [7]-[9]. 

An automatic assessment for blood vessel anomalies of the 
optic fundus initially requires the segmentation of the vessels 
from the background, so that suitable feature extraction and 
processing may be performed. Several methods have been 
developed for vessel segmentation, but visual inspection and 
evaluation by receiver operating characteristic (ROC) analysis 
shows that there is still room for improvement [10], [11]. 
In addition, it is important to have segmentation algorithms 
that do not critically depend on configuring several parameters 
so that untrained community health workers may utilize this 
technology. These limitations of the state-of-the-art algorithms 
have motivated the development of the framework described 
here, which only depends on manually segmented images. 

Many different approaches for automated vessel segmenta- 
tion have been reported. The papers [12]— [18] present vessel 
tracking methods to obtain the vasculature structure, along 
with vessel diameters and branching points. Tracking consists 
of following vessel center lines guided by local information, 
usually trying to find the path which best matches a vessel 
profile model. The use of deformable models also shows 
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promising results in [19]-[22]. In [2], [23], [24], matched 
filters are used to emphasize blood vessels. An improvement 
is obtained in [2] by a region-based threshold probing of 
the matched filter response. Multithreshold probing is directly 
applied to the images in [25]. A non-linear filter that enhances 
vessels by exploiting properties of the vessel profiles is in- 
troduced in [26]. Along this line is the use of mathematical 
morphology filtering in [27], [28], coupled with curvature 
evaluation. In [29], multi-scale curvature and border detection 
are used to drive a region growing algorithm. 

Supervised methods for pixel classification have been shown 
in [1], [30], [31]. In [30], feature vectors are formed by 
gray-scale values from a window centered on the pixel being 
classified. A window of values is also used in [31], but 
the features used are a principal component transformation 
of RGB values and edge strength. In [1], ridge detection 
is used to form line elements and partition the image into 
patches belonging to each line element. Pixel features are 
then generated based on this representation. Many features 
are presented and a feature selection scheme is used to select 
those which provide the best class separability. 

Previously, we have shown promising preliminary results 
using the continuous wavelet transform (CWT) [32], [33] 
and integration of multi-scale information through supervised 
classification [34]. Here we improve on those methods using 
a Bayesian classifier with Gaussian mixture models as class 
likelihoods and evaluate performances with ROC analysis. 
ROC analysis has been used for evaluation of segmentation 
methods in [1], [2], [25] and comparison of some of the cited 
methods in [10], [11]. 

In our approach, each pixel is represented by a feature 
vector including measurements at different scales taken from 
the continuous two-dimensional Morlet wavelet transform. The 
resulting feature space is used to classify each pixel as either 
a vessel or non-vessel pixel. We use a Bayesian classifier with 
class-conditional probability density functions (likelihoods) 
described as Gaussian mixtures, yielding a fast classification, 
while being able to model complex decision surfaces and 
compare its performance with the linear minimum squared 
error classifier. 

Originally devised for suitably analyzing non-stationary and 
inhomogeneous signals, the time-scale analysis took place to 
accomplish unsolvable problems within the Fourier frame- 
work, based on the continuous wavelet transform (CWT). 
The CWT is a powerful and versatile tool that has been 
applied to many different image processing problems, from 
image coding [35] to shape analysis [36]. This success is 
largely due to the fact that wavelets are especially suitable for 
detecting singularities (e.g. edges and other visual features) 
in images [37], extracting instantaneous frequencies [38], and 
performing fractal and multi-fractal analysis. Furthermore, the 
wavelet transform using the Morlet wavelet, also often referred 
to as Gabor wavelet, has played a central role in increasing our 
understanding of visual processing in different contexts from 
feature detection to face tracking [39]. The Morlet wavelet 
is directional and capable of tuning to specific frequencies, 
allowing it to be adjusted for vessel enhancement and noise 
filtering in a single step. These nice characteristics motivate 



the adoption of the Morlet wavelet in our proposed framework. 

This work is organized as follows. The databases used 
for tests are described in Subsection III-AI Subsection III-BI 
presents our segmentation framework based on supervised 
pixel classification. In Subsection IH-CI the feature generation 
process is described, including the 2-D CWT and Morlet 
wavelet. Our use of supervised classification and the classifiers 
tested are presented in Subsection III-DI ROC analysis for 
performance evaluation is described in Subsection III-EI and 
results are presented in SectionED] Discussion and conclusion 
are in Section ITV1 

II. Materials and methods 

A. Materials 

There are different ways of obtaining ocular fundus images, 
such as with non-mydriatic cameras, which do not require 
the dilation of the eyes through drops, or through angiograms 
using fluorescein as a tracer [5]. We have tested our methods 
on angiogram gray-level images and colored non-mydriatic 
images [32], [34]. Here, our methods are tested and evaluated 
on two publicly available databases of non-mydriatic images 
and corresponding manual segmentations: the DRIVE [1] and 
STARE [2] databases. 

The DRIVE database consists of 40 images (7 of which 
present pathology), along with manual segmentations of the 
vessels. The images are captured in digital form from a Canon 
CR5 non-mydriatic 3CCD camera at 45° field of view (FOV). 
The images are of size 768 x 584 pixels, 8 bits per color 
channel and have a FOV of approximately 540 pixels in 
diameter. The images are in compressed JPEG-format, which 
is unfortunate for image processing but is commonly used in 
screening practice. 

The 40 images have been divided into a training and test 
set, each containing 20 images (the training set has 3 images 
with pathology). They have been manually segmented by three 
observers trained by an ophthalmologist. The images in the 
training set were segmented once, while images in the test set 
were segmented twice, resulting in sets A and B. The observers 
of sets A and B produced similar segmentations. In set A, 
12.7% of pixels where marked as vessel, against 12.3% vessel 
for set B. Performance is measured on the test set using the 
segmentations of set A as ground truth. The segmentations 
of set B are tested against those of A, serving as a human 
observer reference for performance comparison. 

The STARE database consists of 20 digitized slides captured 
by a TopCon TRV-50 fundus camera at 35° FOV. The slides 
were digitized to 700 x 605 pixels, 8 bits per color channel. 
The FOV in the images are approximately 650 x 550 pixels in 
diameter. Ten of the images contain pathology. Two observers 
manually segmented all images. The first observer segmented 
10.4% of pixels as vessel, against 14.9% vessels for the second 
observer. The segmentations of the two observers are fairly 
different in that the second observer segmented much more of 
the thinner vessels than the first one. Performance is computed 
with the segmentations of the first observer as ground truth. 
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B. General framework 

The image pixels of a fundus image are viewed as objects 
represented by feature vectors, so that we may apply statistical 
classifiers in order to segment the image. In this case, two 
classes are considered, i.e. vessel x non-vessel pixels. The 
training set for the classifier is derived by manual segmen- 
tations of training images, i.e. pixels segmented by hand are 
labeled as vessel while the remaining pixels are labeled as non- 
vessel. This approach allows us to integrate information from 
wavelet responses at multiple scales in order to distinguish 
pixels from each class. 

C. Pixel features 

When the RGB components of the non-mydriatic images 
are visualized separately, the green channel shows the best 
vessel/background contrast (Fig. |l(a)) , whereas, the red and 
blue channels show low contrast and are very noisy. Therefore, 
the green channel was selected to be processed by the wavelet, 
as well as to compose the feature vector itself, i.e. the green 
channel intensity of each pixel is taken as one of its features. 
For angiograms, the wavelet is applied directly to the gray- 
level values, which are also used to compose the feature 
vectors. 

1 ) Pre-processing: In order to reduce false detection of the 
border of the camera's aperture by the wavelet transform, an 
iterative algorithm has been developed. Our intent is to remove 
the strong contrast between the retinal fundus and the region 
outside the aperture (see Fig. [0- 

The pre-processing algorithm consists of determining the 
pixels outside the aperture that are neighbors to pixels inside 
the aperture and replacing each of their values with the mean 
value of their neighbors inside the aperture. This process is 
repeated and can be seen as artificially increasing the area 
inside the aperture, as shown in Fig. |l(b)| 

Before the application of the wavelet transform to non- 
mydriatic images, we invert the green channel of the image, 
so that the vessels appear brighter than the background. 

2) Wavelet transform features: The notation and definitions 
in this section follow [40]. The real plane 1 x R is denoted 
as M 2 , and the vectors are represented as bold letters, e.g. 
x, b £ R 2 . Let / £ L? be an image represented as a square 
integrable (i.e. finite energy) function defined over M 2 . The 
continuous wavelet transform T^(h, 9, a) is defined as: 

T^(h,6,a)=C- l/2 - a J V*(a-V_ e (x-b))/(x)d 2 x 

where C$, ip, b, 9 and a denote the normalizing constant, 
analyzing wavelet, the displacement vector, the rotation angle 
and the dilation parameter (also known as scale), respectively. 
ip* denotes the complex conjugate of ifi. 

Combining the conditions for both the analyzing wavelet 
and its Fourier transform of being well localized in the time 
and frequency domain plus the requirement of having zero 
mean, one realizes that the wavelet transform provides a local 
filtering at a constant rate — , indicating its great efficiency 
as the frequency increases, i.e. as the scale decreases. This 
property is what makes the wavelet effective for detection and 



analysis of localized properties and singularities [38], such as 
the blood vessels in the present case. 

Among several available analyzing wavelets, for instance, 
the 2-D Mexican hat and the optical wavelet, we chose the 
2-D Morlet wavelet for the purposes of this work, due to 
its directional selectiveness capability of detecting oriented 
features and fine tuning to specific frequencies [38], [40]. 
This latter property is especially important in filtering out 
the background noise of the fundus images. The 2-D Morlet 
wavelet is defined as: 

V>m(x) = exp(jk x)exp ^-i|Ax| 2 ^ 

where j = \f~^\ and A = diag[e~ 1//2 , l],e > 1 is a 2 x 2 
diagonal matrix that defines the anisotropy of the filter, i.e. its 
elongation in any desired direction [38]. The Morlet wavelet 
is actually a complex exponential modulated Gaussian, where 
ko is a vector that defines the frequency of the complex 
exponential. 

We have set the e parameter to 8, making the filter elongated 
and k = [0, 3], i.e. a low frequency complex exponential 
with few significant oscillations, as shown in Fig. [2] These 
two characteristics have been chosen in order to enable the 
transform to present stronger responses for pixels associated 
with the blood vessels. 

For each considered scale value, we are interested in the re- 
sponse with maximum modulus over all possible orientations, 
i.e.: 

M^(b, a) = max |T</,(b, 9, a)\ (1) 

9 

Thus, the Morlet wavelet transform is computed for 9 
spanning from up to 170 degrees at steps of 10 degrees and 
the maximum is taken (this is possible because |T^(b, 9,a)\ = 
|7L,(b, 6 + 180, a)|). The maximum modulus of the wavelet 
transform over all angles for multiple scales are then taken as 
pixel features. M^,(h, a) is shown in Fig. [3] for a = 2 and 
a = 4 pixels. 

3) Feature normalization: Given the dimensional nature of 
the features forming the feature space, one must bear in mind 
that this might give rise to errors in the classification process, 
as the units chosen might affect the distance in the feature 
space. 

A strategy to obtain a new random variable with zero mean 
and unit standard deviations, yielding, in addition, dimension- 
less features, is to apply the normal transformation to the 
feature space. The normal transformation is defined as [36]: 

Vi - m 

Vi = 

where v.- L is the i th feature assumed by each pixel, /i^ is 
the average value of the i th feature and o~i is the associated 
standard deviation. 

We have applied the normal transformation separately to 
each image's feature space, i.e., every image's feature space is 
normalized by its own means and standard deviations, helping 
to compensate for intrinsic variation between images (e.g. 
illumination). 



(a) M^,(b, 2). 



(b) M^,(b,4). 



(a) Surface representation of the real part. 



Fig. 3. Maximum modulus of Morlet wavelet transform over angles, 
M^(b, a) (Eq.Q, for scale values of a = 2 and a = 4 pixels. The remaining 
parameters are fixed at e = 8 and ko = [0, 3]. 




(b) Real part. 



(c) Imaginary part. 



Fig. 2. Different representations for the 2-D Morlet wavelet (»/>m) with 
parameters e = 8 and ko = [0, 3]. 



D. Supervised classification for segmentation 

Supervised classification has been applied to obtain the 
final segmentation, with the pixel classes defined as C\ = 
{vessel pixels} and C2 = {non-vessel pixels}. In order to 
obtain the training set, several fundus image have been man- 
ually segmented, allowing the creation of a labeled training 
set into classes C\ and C2 (see Subsection Ill-Ai . Due to 
the computational cost of training the classifiers and the 
large number of samples, we randomly select a subset of the 
available samples to use for actually training the classifiers. 
We will present results for two different classifiers, described 
below. 

1) Gaussian mixture model Bayesian classifier: We have 
achieved very good results using a Bayesian classifier in 
which each class-conditional probability density function (like- 
lihood) is described as a linear combination of Gaussian 
functions [41], [42]. We will call this the Gaussian mixture 
model (GMM) classifier. 

The Bayes classification rule for a feature vector v can be 
stated in terms of posterior probabilities as 



Decide d if P(C*i|v) > P(C 2 |v); 
otherwise, decide C2 



(2) 



We recall Bayes rule: 
P(Ci\v) 



P (v|q)p(q 

p(v) 



(3) 



where p(v\d) is the class-conditional probability density 
function, also known as likelihood, P{Ci) is the prior proba- 
bility of class Ci, and p(v) is the probability density function 
of v (sometimes called evidence). 

To obtain a decision rule based on estimates from our 
training set, we apply Bayes rule to Eq. |2] obtaining the 
equivalent decision rule: 



Decide C x if p{-v\Ci)P{Ci) 
otherwise, decide C2 



>j>(v|C a )p(C 2 ); 



We estimate P(Ci) as Ni/N, the ratio of class i samples in 
the training set. The class likelihoods are described as linear 
combinations of Gaussian functions: 



p(v\Ci) 



hi 

E 



p(v|j,a)p, 



where fc, is the number of Gaussians modeling likelihood i, 
Pj is the weight of Gaussian j and each p(v\j,Ci) is a d- 
dimensional Gaussian distribution. 

For each class i, we estimate the fej Gaussian parameters 
and weights with the Expectation-Maximization (EM) algo- 
rithm [41]. The EM algorithm is an iterative scheme that 
guarantees a local maximum of the likelihood of the training 
data. 

GMMs represent a halfway between purely nonparametric 
and parametric models, providing a relatively fast classification 
process at the cost of a more expensive training algorithm. 

2) Linear minimum squared error classifier: We have also 
tested the linear minimum squared error classifier [41], [42], 
denoted LMSE. Linear classifiers are defined by a linear 
decision function g in the d-dimensional feature space: 



w 



(4) 
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where v is a feature vector, w is the weight vector and wq 
the threshold. 

The classification rule is to decide C\ if g(v) > and C2 
otherwise. To simplify the formulation, the threshold wq is 
accommodated by defining the extended (d + 1) -dimensional 
vectors v' = [v T ,l] T and w' = [w t ,wq} t , so that g(v) = 



IT 1 

w v . 



The classifier is determined by finding w' that minimizes 
the sum of error squares criterion: 



N 



J(w')=$> 



w') 2 



where N is the total number of training samples, is the 
extended i th training sample, and yi its desired output. 

The criterion measures the sum of squared errors between 
the true output of the classifier (v'^w') and the desired output 
(Ui)- We have arbitrarily set yi = 1 for Vj <E C\ and yi = — 1 
for v, e C 2 . 

Let us define 



V = 



y = 



!)2 
UN 



Minimizing the criterion with respect to w' results in: 

(V T V)w' = V T y => w' = (V T V)- l V T y 

In comparison to the GMM classifier, the LMSE classi- 
fier has a much faster training process, but is restricted in 
the sense that it is linear, while GMMs allow for complex 
decision boundaries. However, as we will show, the results 
obtained using LMSE are comparable to those using GMMs, 
representing a reasonable trade-off. 



E. Experimental evaluation 

The performances are measured using receiver operating 
characteristic (ROC) curves. ROC curves are plots of true 
positive fractions versus false positive fractions for varying 
thresholds on the posterior probabilities. A pair formed by a 
true positive fraction and a false positive fraction is plotted 
on the graph for each threshold value (as explained below), 
producing a curve as in Figs. [6] and The true positive 
fraction is determined by dividing the number of true positives 
by the total number of vessel pixels in the ground truth 
segmentations, while the false positive fraction is the number 
of false positives divided by the total number of non-vessel 
pixels in the ground truth. In our experiments, these fractions 
are calculated over all test images, considering only pixels 
inside the FOV. 

For the GMM classifier, the ROC curve is produced by 
varying the threshold on the posterior pixel probabilities (see 
Eq. [3}, while the LMSE ROC curve is produced varying the 
threshold wq on the projection of the feature vectors on the 
discriminant vector (see Eq. |4}. 



We have tested our methods on the DRIVE and STARE 
databases with the following settings. The pixel features used 
for classification were the inverted green channel and its 
maximum Morlet transform response over angles M^,(b, a) 
(Eq. Q for scales a — 2,3,4,6 pixels (see Subsection III-O . 

For the DRIVE database, the training set was formed by 
pixel samples from the 20 labeled training images. For the 
STARE database, leave-one-out tests where performed, i.e., 
every image is segmented using samples from the other 19 
images for the training set. Due to the large number of pixels, 
in all experiments, one million pixel samples where randomly 
chosen to train the classifiers. Tests were performed with the 
LMSE and GMM classifiers. For the GMM classifier, we vary 
the number k = k\ = k^ of vessel and non-vessel Gaussians 
modeling each class likelihood. 

III. Results 

Illustrative segmentation results for a pair of images from 
each database (produced by the GMM classifier with k = 20), 
along with the manual segmentations, are shown in Figs. 0] 
and|3 

For the DRIVE database, the manual segmentations from 
set A are used as ground truth and the human observer 
performance is measured using the manual segmentations from 
set B, which provide only one true/false positive fraction pair, 
appearing as a point in the ROC graph (Fig.|6). For the STARE 
database, the first observer's manual segmentations are used 
as ground truth, and the second observer's true/false positive 
fraction pair is plotted on the ROC graph (Fig. [7}. The closer 
an ROC curve approaches the top left corner, the better the 
performance of the method. A system that agreed completely 
with the ground truth segmentations would yield an area under 
the ROC curve A z = 1. However, note that the second sets 
of manual segmentations do not produce perfect true/false 
positive fractions, for the manual segmentations evaluated 
disagree on some of the pixels with the manual segmentations 
used as ground truth. Thus, the variance between observers 
can be estimated, helping to set a goal for the method's 
performance. 

The areas under the ROC curves (A z ) are used as a single 
measure of the performance of each method and are shown in 
Table |2 for GMM classifiers of varying k and for the LMSE 
classifier. For comparison with the manual segmentations, we 
also measure the accuracies (fraction of correctly classified 
pixels) of the automatic and manual segmentations. Note that 
the accuracy and A z values for the GMM classifier increase 
with k. The ROC curves for the DRIVE and STARE databases 
produced using the GMM classifier with k = 20, as well as 
performances for human observers, are shown in Figs. |6] and [7J 

We note that the EM training process for the GMMs is 
computationally more expensive as k increases, while the 
classification phase is fast. On the other hand, LMSE is very 
fast for both training and classification, but produces poorer 
results, as seen in Table ||] 

IV. Discussion and conclusion 

The Morlet transform shows itself efficient in enhancing 
vessel contrast, while filtering out noise. Information from 



(a) Posterior probabilities. 



(b) Segmentation. 



(c) Set A. 



(d) Set B. 




(e) Posterior probabilities. (f) Segmentation. (g) Set A. (h) Set B. 

Fig. 4. Results produced b y the GMM classifier with k = 20 and manual segmentations for two images from the DRIVE database. The top row results are 
for the image shown in Fig. |l(a)| 




(a) Posterior probabilities. (b) Segmentation. (c) First observer. (d) Second observer. 




(e) Posterior probabilities. (f) Segmentation. (g) First observer. (h) Second observer. 

Fig. 5. Results produced by the GMM classifier with k = 20 and manual segmentations for two images from the STARE database. The top row images 
originate from a pathological case, while the bottom ones originate from a normal case. 
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2 0.6 



GMM, k = 20 

□ Set B 



false positive fraction 



Fig. 6. ROC curve for classification on the DRIVE database using the 
GMM classifier with k = 20. The point marked as □ corresponds to set B, 
the second set of manual segmentations. The method has A z = 0.9598. 




GMM,k = 20 

□ 2nd oOserver 



0.2 0.3 0.4 0.5 0.6 0.7 

false positive fraction 



0.8 0.9 



Fig. 7. ROC curve for classification on the STARE database using the 
GMM classifier with k = 20. The point marked as □ corresponds to the 
second observer's manual segmentations. The method has A z = 0.9651. 



TABLE I 

Results for different classification methods and human 
observer. a z indicates the area under the roc curve, while 
the accuracy is the fraction of pixels correctly classified. 



Classification Method 


Database 


DRIVE 


STARE 


A z 


Accuracy 


A z 


Accuracy 


LMSE 


0.9520 


0.9280 


0.9584 


0.9362 


GMM, k = 1 


0.9250 


0.9217 


0.9394 


0.9239 


GMM, k = 5 


0.9537 


0.9431 


0.9609 


0.9430 


GMM, k = 10 


0.9570 


0.9454 


0.9627 


0.9450 


GMM, k = 15 


0.9588 


0.9459 


0.9648 


0.9470 


GMM, k = 20 


0.9598 


0.9467 


0.9651 


0.9474 


2nd. observer 




0.9473 




0.9349 



Morlet transforms at different scales, which allows the seg- 
mentation of vessels of different diameters, are integrated 
through the use of the statistical classifiers presented. The 
LMSE classifier shows a reasonable performance with a fast 
classification and training phase, while the GMM classifier has 
a computationally demanding training phase, but guarantees a 
fast classification phase and better performance. 

The classification framework demands the use of manual 
labelings, but allows the methods to be trained for different 
types of images (provided the corresponding manual segmen- 
tations are available), possibly adjusted to specific camera or 
lighting conditions and are otherwise automatic, i.e., adjust- 
ment of parameters or user interaction is not necessary. We are 
studying the use of training sets composed of a small portion 
of the image to be segmented. Using this approach, a semi- 
automated fundus segmentation software may be developed, 
in which the operator only has to draw a small portion of the 
vessels over the input image or simply click on several pixels 
associated with the vessels. The remaining image would then 
be segmented based on the partial training set. This approach 
is interesting since it requires a small effort from the operator, 
which is compensated by the fact that image peculiarities are 
directly incorporated by the classifier. 

It is curious to note that, on the STARE database, the 
accuracy of the method is higher than that of the second 
observer (Table [Q. The second observer's manual segmenta- 
tions contain much more of the thinnest vessels than the first 
observer (lowering their accuracy), while the method, trained 
by the first observer, is able to segment the vessels at a similar 
rate. However, the ROC graph (Fig. still reflects the higher 
precision of the second observer, due to some difficulties found 
by the method, as discussed below. 

It is possible to use only the skeleton of the segmentations 
for the extraction of features from the vasculature. Depending 
on the application, different evaluation methods become more 
appropriate [43]. For example, the evaluation of the skeleton 
would not take into account the width of the vessels, but 
could measure other qualities such as the presence of gaps 
and detection of branching points. Another interesting form 
of evaluation would be directly through an application, such 
as in detection of neovascularization by means of analysis and 
classification of the vessel structure [33]. A major difficulty in 
evaluating the results is the establishment of a reliable ground 
truth [44]. Human observers are subjective and prone to errors, 
resulting in large variability between observations. Thus, it 
is desirable that multiple human-generated segmentations be 
combined to establish a ground truth, which was not the case 
in the analysis presented. 

Though very good ROC results are presented, visual inspec- 
tion shows some typical difficulties of the method that must be 
solved by future work. The major errors are in false detection 
of noise and other artifacts. False detection occurs in some 
images for the border of the optic disc, haemorrhages and 
other types of pathologies that present strong contrast. Also, 
the method did not perform well for very large variations in 
lighting throughout an image, but this occurred for only one 
image out of the 40 tested from both databases. This could 
possibly be solved by including intra-image normalization 
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in the pre-processing phase [45]. Another difficulty is the 
inability to capture some of the thinnest vessels that are barely 
perceived by the human observers. 

Another drawback of our approach is that it only takes 
into account information local to each pixel through image 
filters, ignoring useful information from shapes and structures 
present in the image. We intend to work on methods addressing 
this drawback in the near future. The results can be slightly 
improved through a post-processing of the segmentations for 
removal of noise and inclusion of missing vessel pixels as 
in [34]. An intermediate result of our method is the intensity 
image of posterior probabilities, which could possibly benefit 
from a threshold probing as in [2] or region growing schemes. 

Automated segmentation of non-mydriatic images provides 
the basis for automated assessment by community health 
workers. Skeletonized images of the vessel pattern of the 
ocular fundus can be analyzed mathematically using nonlinear 
methods such as global fractal [33] and local fractal [7] anal- 
ysis based on the wavelet transform thus providing a numeric 
indicator of the extent of neovascularization. Our ongoing 
work aims at applying the shape analysis and classification 
strategies described in [33] to the segmented vessels produced 
by method described in this work. 
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